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Abstract 

We revise the basic concepts beneatli the idea of superparamagnetism 
and the suitability of Monte Carlo (MC) simulations to study superparam- 
agnetic (SPM) properties. Starting with the description of the characteris- 
tic features of the single-domain SPM entities, their general magnetic-field 
and temperature-dependent magnetic properties are discussed. Then, the 
use of a MC technique for studying SPM properties is presented, starting 
with a general approach to MC methods and introducing the Metropolis 
algorithm as an adequate tool for reproducing SPM features. Special at- 
tention is paid to the role played by the computational time MC steps on 
the simulations. 



1 Introduction 

Nanosized magnetic materials exhibit a rich variety of magnetic phenomena 
in comparison with the bulk counterparts, what gives place to a novel range 
of appUcations of great importance to improve daily human activities as high- 
density magnetic recording storage or biomedical applications ITj. The origin 
of these special magnetic properties is found on the reduced dimensionality [2]: 
when the size of the material reaches the order of nanometers the influence 
of the surface atoms becomes very comparable (or even higher) than the bulk 
contribution, the defects due to the broken symmetry may be of significant 
importance, and other physical effects may also become very relevant when 
the size reaches the order of characteristic length scales of the material (as for 
example the domain size) . The properties observed in such reduced dimensions 
are strongly sensitive to small size, shape, and composition dependence, what 
defines the different magnetic structures (nanoparticles, nanoparticle arrays. 
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nanowires, thin films, etc) as forming specific research fields with differentiated 
proper characteristics [3]. 

In the first part of this article we revise one of the most remarkable magnetic 
properties that arises in these reduced dimensions, the so-called superparamag- 
netism. Superparamagnetic (SPM) phenomena is the paramagnetic-like tem- 
perature dependence that occurs in single-domain magnetic entities above a 
characteristic threshold named blocking temperature, and its special features are 
determined by the complex interplay between the magnetic parameters ruling 
the system (magnetic moment, anisotropy, applied field, etc). Understanding 
the SPM properties of nanosized systems is of primordial importance both for 
the basic theoretical knowledge \4\ and for specific-designed applications (as for 
example the increase of the storage information capacity of hard drives [5] , or 
the development of well-controlled biomedical applications [6]). Because of this 
a big effort has been devoted to the study of SPM systems in the last years, 
aimed to understand its underlying physical mechanisms. 

However, the investigation of SPM properties is a complex task due to its 
strong dependence on several uncontrolled parameters, which mask the physical 
origin of the magnetic behaviour and hence makes very difficult to achieve a pre- 
cise characterization. The high parameter-dispersion degree found in real sys- 
tems, arising from the large dispersion of parameters (particle size, anisotropy, 
shape) and uncontrolled interparticle interactions, results in a complex physi- 
cal problem non-solvable by analytical methods. To investigate such complex 
scenario it is very common the use of computational techniques, which allow a 
precise control of the physical parameters governing the system: by means of 
a computational technique it is possible to set up ideal systems (e.g. monodis- 
perse size and/or anisotropy) specially designed to ease the comprehension of 
the physical mechanisms governing the system. In the second part of this work 
we introduce the basic characteristics of a a Monte Carlo (MC) method based 
on the Metropolis algorithm to undertake the study of SPM properties. 

2 Superparamagnetism 

The term superparamagnetism refers to the magnetic phenomena observed in 
fine magnetic particle systems exhibiting close similarities to atomic paramag- 
netism. Basically, single-domain magnetic nanoparticles can be characterized 
by their large total magnetic supermoment, which exhibits a paramagnetic-like 
decay of the magnetization with temperature above a characteristic threshold 
named blocking temperature, Tb- This particular temperature, as difference to 
the Curie temperature Tq, is extremely dependent on the experimental observa- 
tional time-window and this characteristic gives place to a complex theoretical 
frame with especial experimental features. In what follows we briefly introduce 
the conditions for the existence of superparamagnetism and its basic charac- 
teristics. For more information in this topic see for example the reviews by D. 
KechrakosQ, O. Petracic[8], M. Knobel et al.[9^, J.L. Dormann et al.[4]. 
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2.1 Single-domain particles 



SPM phenomena is observed upon reduction of the size of ferromagnetic (FM) 
materials|10| down to the single-domain range: in a FM magnetic material, mul- 
tiple magnetic domains exist as a result of the balance between the exchange 
interaction energy favouring the parallel alignment of neighbouring atomic mo- 
ments (thereby forming magnetic domains), and the magnetostatic interaction 
energy forcing their breaking into smaller domains with tendency to antiparal- 
lel orientation. The domain size is determined by the relative counterbalance 
between both energies. With decreasing size of the magnetic system, there is 
a critical value (rc for the radius of a spherical particle) below which the mag- 
netostatic energy no longer allows for the breaking of the system into smaller 
domains and so the system is composed of a single domain, as illustrated in 
Figure [T] Typical values are Vc ~ 15 nm for Fe and « 35 nm for Co [1]. 



Figure 1: Scheme illustrating the transition from the multi-domain configuration 
to the single-domain one upon size reduction. 

Assuming coherent rotation of the domain atomic moments, the particle is 
therefore characterized by its total magnetic supermoment, flp, resulting from 
the total magnetization of the particle. In first approximation, considering 
uniform magnetization (so neglecting surface effects) it can be described as 
proportional to the particle volume V and saturation magnetization Ms as 



As mentioned above, the SPM response in magnetic nanoparticles is ob- 
served above the so-called blocking temperature Tb, a proper feature of SPM 
systems that differentiates them from atomic paramagnetism. The origin of Tb 
relays on the magnetic anisotropy present in the nanoparticles (in opposition 
to atomic moments) due to their finite size, which tends to orientate the par- 
ticle supermoment along some preferential direction. The magnetic anisotropy 
energy Ea found in a magnetic nanoparticle can have different origins (crystal, 
shape, surface, etc) giving place to very complex scenarios, and so for the sake 
of simplicity we have focused on the simplest uniaxial anisotropy casepj]. So, 
from now on we consider the different anisotropy contributions as comprised in 
an effective uniaxial anisotropy term, Kce, as illustrated in Figure [Sf a). 

If we consider the magnetic anisotropy to be proportional to the particle vol- 
ume as A'cff — KVh, with K the effective uniaxial anisotropy constant (per unit 




Ia^pI = MsV 



(1) 
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volume) and h the unitary vector describing the easy-magnetization anisotropy 
direction, then the energy term for the z-particle can be written as 

= -KV, f ' = -KV^ cos2 e (2) 

being 9 the angle between the magnetic supermoment of the particle and the 
easy anisotropy axis. The moment of the particle has therefore two preferred 
orientation, equally probable, along the easy-magnetization anisotropy axis di- 
rection. Both directions are separated by an energy barrier Eb of height KV . 
The energy spectra corresponding to this uniaxial anisotropy energy is illus- 
trated in Figure [2jb). 
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Figure 2: (a) Schematic drawing of the uniaxial magnetic anisotropy K and mag- 
netic supermoment of a single-domain nanoparticle, and (b) the corresponding 
uniaxial anisotropy wells. 

The system we have constructed up to now is that of homogeneous magnetic 
nanoparticles characterized by their size V ^ saturation magnetization Ms, and 
uniaxial magnetic anisotropy energy K . This very simple scenario describes rea- 
sonably well many experimental situations, and so from now on we focus on the 
magnetic properties of a system of such particles as a function of temperature 
(T) and magnetic field {H). Real systems are usually characterized by random- 
ness in their spatial distribution and in the easy-axes orientation that strongly 
influence the properties of the system as determined by interparticle- interactions 
and applied magnetic field. These analysis are however quite complex and so 
for the sake of simplicity we consider for the moment a non-interacting parti- 
cle system (so spatial distribution concerns are avoided) with parallel aligned 
anisotropy axes, so that the particles are equivalent to each other and the system 
can be studied under a single-particle perspective. 
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2.2 Thermal relaxation and blocking temperature (Tb) 

The influence of the temperature on the magnetic properties of the particles 
can be easily figured out in the high- and low-T limit cases. At very high 
T, the thermal energy is much larger than the anisotropy barrier (fc^T >> 
Eb, with fcs the Boltzmann constant) and so the magnetic anisotropy plays a 
negligible influence on the orientation of the magnetic moments of the particles. 
In this case a paramagnetic-like dependence of the magnetization with T is 
expected and the particles are in the superparamagnetic state (SPM state). On 
the contrary, at very low T the particle moment remains confined along the 
anisotropy direction (local energy minimum) because the thermal energy being 
unable to switch its orientation out of that minimum {ksT << Eb)- When this 
happens the particles are said to be in the blocked state (BL state). 

Thermal energy promotes the fluctuation of the magnetic moments, its in- 
fluence varying between the rapid motion at very high temperatures and the 
practically steady state at very low T. Hence, to understand how thermal en- 
ergy influences the magnetic behavior of the particles it is necessary to under- 
stand the dynamics of the particle moments as a function of T. The influence of 
thermal fluctuation on the orientation of the particles' supermoments was first 
described by Neel [[12 , who proposed that the thermal fluctuations could pro- 
mote the jumping of the magnetic moment of the particles from one anisotropy 
well to the other, introducing the average time r for thermal activation (often 
called relaxation time) over the anisotropy barrier to follow an Arrhenius law 

T = Toe ''bT (3) 

where tq directly depends on the material parameters {K, Ms, etc) and is of 
the order of 10"^"'^ — 10~®s. Under this description, it therefore points out that 
the measuring time will be a key-point on determining the magnetic state of 
the system: if the measuring time is large in comparison with the characteristic 
reversal time of the particles, Tm >> t, then the particle moment will fluctuate 
rapidly from one well to the other in a paramagnetic-like fashion, i.e. the particle 
will be in the SPM state. However, if Tm << t, during the measuring time the 
particle moment will remain blocked along one anisotropy well, i.e. the particle 
will be in the BL state. Macroscopically, the SPM state is completely reversible 
upon temperature and field variations, whereas the BL one is characterized 
by its hysteretic behavior, proper of ferromagnets. The limit between both 
states is found at = t, and serves for the definition of the so-called blocking 
temperature Tb , as illustrated in Figure [3] and obtained from Eq. ([3]) 

KV 

Tb = (4) 

Kb m (Tm/To) 

As Eq.(|3]) shows, Tb depends both on the intrinsic particle parameters and 
on the external ones as the measuring time. Therefore, by varying the exter- 
nal influences of the systems (temperature, measuring time, applied field) we 
may tune its response. From the trends obtained we can extract information 
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about the characteristics of the system, and so information for the design of 
technological applications. 
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Figure 3: Relaxation time r vs T, and indication of Tb for a certain Tm- 

Eq.Q highlights that the value of t„i is decisive in determining Tb and so the 
magnetic response of the system in a given time-scale (for example, for informa- 
tion storage-related applications, very large time scales have to be considered; 
however, for magnetic recording speed very short ones must be achieved). Since 
Tm is determined by the experimental technique, its value has to be chosen with 
respect to the information and uses concerning our purposes. In this work we 
have focused on quasi static process, and so associate Tb with the one obtained 
in dc-thermomagnetization measurements, in which the measuring time is very 
large large (r^ ~ 100s) in comparison with the characteristic time tq of the par- 
ticles. Other measurement techniques involving much shorter measuring times 
are associated to dynamic measurements, not considered in this work [14]. We 
have mainly focused on data obtained following the standard zero field cooling 
(ZFC) and field cooling (FC) protocols, in which the system is perturbed un- 
der a low magnetic field for measuring the evolution of its magnetization with 
temperature. Previous to describing such processes and its characteristics, we 
shall analyze the effect of the magnetic field on the magnetic properties of the 
single-domain nanoparticles. 

2.3 Field dependence 

When an external magnetic field is applied over the nanoparticles, it tries to 
align their magnetic moments along its direction. Therefore, except if applied 
perpendicularly to the easy anisotropy axis, it will favour the occupancy of one 
of the anisotropy energy wells over the other. The orientation of the magnetic 
moment of a particle i is then governed by the competition among its uniaxial 
{Ea) and Zeeman {Ez) energies 
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E^=E^^+E^^^ = -KV,(^^j -fI,-H (5) 

The influence of the external field in the orientation of the magnetic super- 
moments, assuming coherent rotation of the atomic magnetic moments and the 
field applied at a certain angle 9o with respect to the easy anisotropy axis, is 
known as the Stoner- Wohlfart model [TS] after the authors who first considered 
and solved this problem. They ignored thermal effects and so could solve it from 
minimal energy arguments, not taking into account time-dependence as related 
to temperature. In this introduction section, however, we are mainly interested 
in giving an overall description of the SPM features, as how the magnetic field 
will affect the orientation of the magnetic moments with temperature in relation 
to the anisotropy energy barrier and its influence on Tb- Therefore, we do not 
discuss now the orientation of the field at different angles and focus, for the 
sake of simplicity, on the simple case of the field applied parallel to the easy 
anisotropy axis (H\\hi). Note that in this context of non-interacting particles, 
applying the magnetic field at a certain angle with respect to the easy axis is 
equivalent to consider only its projection along the axis. Since we are under 
the assumption of non-interacting and equivalent particles, we can apply single- 
particle considerations and simplify the i subindex in Eq.([5]), which taking into 
account Eq.([T]), reads 

E = ~KV cos^ 9 - MsVH cos 9 (6) 

Eq.([6]) has two local energy minima at = 0, tt with values Emin = —KV ± 
MsVH, and a maximum a. 9 = tt/2 with value E,nax = KV{HMs/2K)^. The 
9 = value stands for the parallel orientation of the particle moment with 
respect to the magnetic field (tt)i whereas the 9 = n value stands for the 
antiparallel one (ti)- This difference in the energy wells described by Figure [2] 
corresponds therefore to different energy barriers depending on the orientation 
of the particle moment with respect to the applied field, which we shall call E^^ 
and e]J for the antiparallel and parallel cases, respectively. Introducing the 
anisotropy field of the particles as 

we calculate these energy barriers as the difference between the minima and 
maximum energies, obtaining 



H 



2 



El,^ = KV{l- — ] (8) 



and 



H 



E^s^ =KV{1 + —] (9) 
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The difference in the height of the energy barriers indicates also a change in 
the characteristic relaxation time of the particles, since it depends now on the 
relative orientation of the magnetic dipoles with respect to the field: particles 
antiparallel-oriented with respect to the field have a smaller energy barrier in 
comparison to the solely-anisotropy one and so a smaller thermal energy is 
enough to overcome it, whereas the parallel-oriented particles are now confined 
into a deeper anisotropy well and so a higher thermal energy is necessary to 
promote the jumping of its magnetic moments. This influence of the magnetic 
field on the energy wells of Figure [5] is shown in Figure g] (left panel) , as well 
as its implications into the relaxation time of the particles (right panel) . 




Figure 4: Anisotropy energy wells (left panel) and relaxation time (right panel) 
of the particles as influenced by the magnetic field. 

Figure d] illustrates the importance of the strength of the applied magnetic 
field on determining the magnetic properties of SPM systems, and also serves 
as a definition of small field as referred to magnetic nanoparticles in comparison 
with their anisotropy field: the ratio H/Ha must be as smaller as possible so 
that the system keeps as closer as possible to the ideal SPM conditions. 

2.4 Thermomagnetic measurements 

Once we have gone through how a magnetic field H influences the properties 
of the particles we can undertake the description of the ZFC and FC measure- 
ments. In both processes the temperature evolution of the total magnetization 
of the system is recorded following different thermomagnetic histories, and it 
is this different history what highlights reversibility (no-hysteresis) and irre- 
versibility (hysteresis) for differentiating the anhysteretic SPM state from the 
hysteretic BL state[Tn]. In a ZFC process, the system is flrst cooled down in 
zero fleld until a very low T is reached, and afterwards a small fleld is applied 
and the magnetization is recorded while heating the sample up. The FC curve 
is obtained by measuring the magnetization of the sample while cooling under 
low magnetic fleld (same field strength for both ZFC and FC processes, and low 
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to minimize the disturbance of the system). Typical ZFC/FC curves are shown 
in Figure [5] 
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Figure 5: Typical ZFC and FC magnetization curves vs temperature. Vertical 
dotted line stands for the maximum of the ZFC curve, usually associated to Tg. 



The ZFC and FC curves shown in Figure [S] display the usual features found 
in SPM systems: i) both curves coincide at high temperatures in a PM-like de- 
pendence; ii) with decreasing T both curves grow until a certain temperature 
range is reached where the curves start to diverge, the FC curve still growing 
although at a lower rate while the ZFC exhibits a maximum and decreases be- 
low it. This maximum in the ZFC curve is generally associated in the literature 
to Tb [7], as indicated in the figure, since such maximum roughly differenti- 
ates two main temperature regimes: a high-temperature one where both curves 
essentially coincide and exhibit a PM-like temperature dependence, from a low- 
temperature regime where both curves clearly diverge. However, a detailed view 
of the curves reveals that a 1 /T PM-like decrease right above Tb is not observed 
in the ZFC curve, and a slight difference between the ZFC and FC one is per- 
ceived. These features indicate that a true SPM behavior is not exhibited right 
above the maximum, but only at higher temperatures the ZFC curve perfectly 
overlaps the FC one and exhibits a well-defined PM-like temperature depen- 
dence. The reason why the ideal SPM behavior is not observed right at T > Tb 
is found on the inverse of the relaxation time r, which gives the probability of 
the particle to overcome the anisotropy energy barrier along the temperature 

p{T) = T-^e'^ (10) 

Eq. pmi considering the different energy barriers E'^ and and so dif- 
ferent relaxation times (see Figure 2]), indicate that SPM behavior will be 
only observed (ideally) above T^^, which above the maximum of the curve. 
Therefore, the features displayed in Figure 2] concerning the different height of 
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the anisotropy wells as influenced by the external field, must be taken into ac- 
count too when analyzing the physical trend followed in the ZFC magnetization 
curve. The initial state of the ZFC process starts with no net magnetization 
after cooling in zero field. If naming parallel particles those with Ez < 0, and 
antiparallel particles those with Ez > 0, then when the field is applied par- 
allel particles will rapidly align with the field, while the antiparallel ones will 
progressively overcome the energy barrier with the increasing thermal energy 
and also align with the field. This process leads to a continuous increase of the 
magnetization, as illustrated in Figure [SJ until the thermal energy overcomes 
E^^, and so no longer reversal magnetization takes place. This is not so simple 
however, since the thermal energy does not induce the alignment of more par- 
ticles, but counterbalances their orientation, inducing thermal fluctuations. As 
the thermal energy is now comparable to the energy of the deeper anisotropy 
well, inverse reversal mechanism overcoming has a higher appreciable proba- 
bility to occur. This frame stands so for the expected further decrease in the 
magnetization, although not yet accomplished full PM-like behavior. In fact, 
pure PM-like behavior would only be observed for fc^T >> Eb- It is important 
to emphasize these aspects when dealing with the magnetic properties of SPM 
systems, remarking that Tb defines only a characteristic temperature value, of 
enormous interest for characterizing the system but not equivalent to a phase 
transition [13 . Because of the above reasons, some authors do not associate 
directly the maximum in the curves to Tb and prefer a different definition [S]. 
Although we understand such discrepancy and share the necessity of finding a 
more precise formalism for relating the shape of the ZFC curves to the particles' 
characteristics, we associate in this work the maximum of the ZFC curve with 
Tb because it is the usually preferred when dealing with experimental data. 

It is important to remind here that all the properties described up to now 
correspond to the simplest single-particle scenario, where all the particles have 
been treated as equivalent to each other because of been non-interacting and 
with parallel-aligned easy anisotropy axes, as illustrated in Figure [SI 

However, real systems are usually characterized by a random distribution 
of the anisotropy easy axes, what would be equivalent to a distribution of the 
heights of the energy barriers along the field direction, as schematized in Fig- 
ure [7| giving place to a extraordinarily complex scenario. Moreover, this situa- 
tion becomes even more complex in real systems, for one of the main problems 
concerning the magnetic properties of SPM systems is the role played by in- 
terparticle dipole-dipole interactions, which are long-range and anisotropic, re- 
sulting so in a complex interplay with the anisotropy barriers. Low-interacting 
conditions can be described by mean-field approximations in which the single- 
particle barriers are slightly modified by the interaction energy. However, strongly- 
interacting conditions ruled by collective effects cannot be accounted by that 
approach, and so the use of computational techniques becomes an indispens- 
able tool for treating systems with so many freedom degrees. Computational 
techniques allow us to deal with perfectly controlled systems and the exact treat- 
ment of the interactions among particles. There are two main computational 
approaches for dealing with the magnetic properties of interacting nanoparticle 
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Figure 6: Schematic drawing of the ideal simplest model of non-interacting and 
parallel aligned easy axes along the applied field. 

systems, namely the Monte Carlo (MC) and Langevin Dynamics (LD) methods. 
Both methods are complementary for the study of a nanoparticle system: MC 
simulations are very adequate to treat long-time (static) magnetic properties in 
complex interacting systems [18] but do not have associated a physical time; 
LD methods, on the contrary, are very precise for modeling the dynamics of 
the magnetic moments but cannot treat static processes.}^ As mentioned 
above we have focused our description of the SPM phenomena by means of its 
features throughout static measurements; next we introduce the MC technique 
in the context of studying complex SPM systems as the ones described here. 

Finally, it is worthy to recall again the several simplifications assumed in 
this introduction to superparamagnetism, where we have considered very simple 
and ideal characteristics for the particles. In real systems there are always 
several dispersive -often uncontrolled- characteristics (inhomogeneities in the 
particles' composition; temperature-dependent K and/or Ms; size/anisotropy 
distribution; etc) , characteristics all that confer additional uncontrolled degrees 
of freedom to the already very complex system, and so make very hard to 
interpret the magnetic measurements in order to characterize their properties. 
That is the reason why we had to focus on a very simple scenario which, on the 
other side, it is already complex enough to make absolutely necessary the use 
of computational techniques for its study. 

3 Monte Carlo method 

Monte Carlo (MC) methods are a type of numerical simulation techniques based 
on the generation of random numbers [23] . MC methods are utilized to solve 
complex problems with large freedom degrees: the features of a particular prob- 
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Figure 7: Schematic drawing of the usual random easy axes distribution, and 
the modulation of the energy barrier height along the field direction. 

lem are represented by probabilities, and the MC technique consists on generat- 
ing large amounts of random numbers and counting the fraction of them obeying 
some conditions. The way of counting and the conditions imposed define the 
numerical algorithm. A simple example to illustrate the functioning of a MC 
method is the calculus of tt from the area of a circle. If placing a circle of radius 
r into a square of side 2r and randomly generating A'^-points into the square, the 
fraction of them laying inside the circle {ricircie) will be equal to the proportion 
between areas, and so it is easy to obtain tt — 4 "°'^^*='^ . The MC calculation of 
the area will be more precise the larger the amount of points (events) generated, 
as illustrated in Figure [H 




Figure 8: Use of random events to calculate the value of tt from the area of a 
circle, and the importance of the amount of events in the precision in the result. 

The two drawings in Figure [8] stand for two examples of random generation 
of events (N=100 and N=1000) into the square with the circle held inside. The 
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graph shows the fraction between points laying in each geometrical figure as 
a function of the amount of random events generated, pointing out that the 
precision to determine tt from the random points areal ratio grows rapidly with 
N. 

3.1 Metropolis algorithm 

We aim in this work to describe how to use a MC method to simulate the 
magnetic properties of a magnetic nanoparticle system with a large amount of 
freedom degrees as described in the previous Section [2l Specifically, we want 
to know how the orientation of the magnetic moments of the particles evolves 
as a function of different parameters (temperature, applied magnetic field, etc), 
i.e. how they behave as a function of the different energies involved. To simulate 
such processes we impose to the system some known conditions and determine 
its configuration from a random generation of events as evaluated under those 
conditions. The processes we want to simulate are essentially quasi static, and 
so the conditions ruling the system can be based on minimum energy arguments 
on the following manner: i) the energy of the system under certain conditions 
is evaluated, ii) an external parameter is varied and the energy is reevaluated, 
and iii) the difference in energy is used to construct a probability function, and 
the change of configuration of the system is accepted or not from the random 
generation of events applied to such probability under a given algorithm. 

It turns out that this problem is much more complex than calculating tt 
from the ratio between the areas of the circle and square as described above, 
where the randomly generated events are equally probable and so the algorithm 
for solving the problem is just to count 1 if the points lay inside the circle 
and otherwise. If we apply the same procedure to simulate the orientation 
of the magnetic moments with temperature, i.e. if we randomly generate new 
possible orientations and evaluate its feasibility to occur, we will find that most 
of the trials are highly improbable and so rejected, and only those with energy 
comparable to the previous state will have some chance of being accepted. For 
example, for simulating the new possible orientation of a particle' magnetic 
moment, initially at an angle Oi (Eq.®), the new trial configuration Of can be 
chosen totally at random, unrelated to 9^, or by considering a slight variation 
after the actual configuration, so that Of = Oi + SO, with 60 small. In the first 
case many trials will be very unfavorable and so rejected, while for the latter 
a higher acceptance ratio is expected. It becomes therefore crucial, in order to 
avoid the wasting of computational time and for optimizing the simulation, to 
be able to select the new trial configurations among the most likely probable 
paths. This can be done if considering a Markov chain of events (configuration 
of one state depends only on the previous one), with the trial state being close 
in energy to the current one. 

The key-points for treating the present problem are, therefore, i) the selec- 
tion of the trial configurations in an efficient way, and ii) the choosing of an 
adequate implementation of the change from the initial state with energy Ei to 
the trial next state with energy Ef. Assuming classical Boltzmann distribution. 
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the probability of a magnetic moment to have energy E at a temperature T is 
proportional to exp{—E/kBT), i.e. p{E) oc e~^/'^«^. Consequently, if consider- 
ing the orientation of the magnetic moments to be markovian the evolution from 
state Ei to state Ef will be proportional to the rate between final and initial 
states probabilities, n^f = p{Ef)/p{Ei) = e^^'^/''^'^, with AE = {Ef - E,). 
This way of choosing the possible next configuration of the system as being en- 
ergetically close to the actual one is named importance sampling, and is based 
on the detailed balance reversibility condition. This approach works very well 
for describing quasi-static thermodynamic processes, as intended in this work, 
although much care has to be taken if dealing with dynamic properties. For 
further details about this topic, see for example: O. Iglesias Doctoral Thesis 
[24], Chapter 5. 

The motion of the magnetic moment of a nanoparticle from the initial state 
with energy Ei to the final state with energy Ef \s often described by means of 
the Metropolis algorithm [2 5) : if Ai? < (the new configuration is more stable 
energetically), the move to the new state is accepted and its energy changes 
to Ef, whereas if AE > (new configuration more unstable), although small, 
the move has still some probability g-^E/kBT occur[26]. To compute this 
probability a random number n with value between and 1 is generated and, 
if n < ri^f the new conformation of the system is accepted and so it has now 
energy Ef, while if n < r,;^/ it is rejected and the energy remains still Ei. 
The Metropolis algorithm for the probability of a AE configuration change is 
expressed as 

(11) 

The Metropolis MC method can be used to simulate the evolution of the 
magnetization of a system of magnetic nanoparticles as a function of different 
parameters. We describe here a MC method based in this procedure (see for 
example Refs. [32l |33l |38l [37l |35l |36l |34]): the simulated system consists on 
an assembly of N-particles contained inside a unit cubic cell (side L), which is 
replicated by using periodic boundary conditions in order to resemble a large 
and homogeneous system. The simulations are always done in two parts: in the 
first one the spatial distribution of the particles is set, and in the second part 
the particles are characterized by their attributes (size, anisotropy, magnetic 
moment) and the evolution of the system is simulated as a function of the 
desired parameters. The positions set in the first part are kept fixed during the 
simulation of the magnetic properties. Next we give a brief description of the 
generalities of the MC method. 

3.2 Spatial arrangement 

The procedure used to achieve the spatial distribution of the particles varies 
depending on the type of system we want to simulate: for ordered structures 
(chains of particles, crystalline structures, etc) the particles are directly placed 
into the unit cubic cell under the desired structure, whereas for setting the 
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spatial arrangement of liquid-like systems (e.g. a frozen fcrrofluid) a relaxation 
algorithm has to be used. 

In Figure [S] we show some chain-like structures obtained by directly plac- 
ing the particles into desired regular positions. They resemble one-dimensional 
columnar parallel chains of magnetic nanoparticles under different spatial ar- 
rangement (square, hexagonal) and different lengths along the X-axis. The 
study of such type of structures is at the center of much research nowadays for 
the basic study of the competition between the enhanced anisotropy and mag- 
netostatic interactions HOI EI] • The magnetic properties of such chain-like 
systems exhibit a good analogy to the behavior of ferromagnetic nanowires [H] . 




Figure 9: Scheme of parallel-aligned chain-like structures hexagonally and 
squared distributed, and different lengths. 

For simulating disordered systems as ferrofluids with liquid-like structure [43j 
the positions of the particles are not directly generated and so we use a Lennard- 
Jones pair potential (ulj) to distribute the particles. During the simulation the 
particles can move freely in space, but their trial positions are markovian-linked 
to the actual one and so the liquid-like structure is more quickly obtained. 
An example of liquid-like structure is shown in FIGURE [TUl together with the 
corresponding correlation function. 

For treating the long-range dipolar interactions the Ewald summation is 
used as in Ref. [18], considering for simulating long and homogeneous systems 
periodic boundary conditions over repetition of the unit cubic cell. The positions 
of the particles and their relative interdistances are calculated and storage for 
the next part of the simulations, in which remain constant throughout. 
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Figure 10: Liquid-like distribution of N=1000 particles into a unit cubic cell and 
correlation function g2{r). Both the correlation function and Lennard- Jones pair 
potential are also indicated. 



3.3 Superparamagnetic (SPM) properties 

Once the spatial distribution of the particles is achieved the next step is to 
characterize them with their main physical properties, namely volume, magnetic 
moment and magnetic anisotropy. Following the model described in Section [2l 
magnetic anisotropy is considered of uniaxial type, and both magnetic moment 
and anisotropy are assumed to be proportional to the particle volume, and so 
the important parameters to characterize the particles are their volume and the 
orientation of the magnetic moment and easy-anisotropy axis. 

The volume is taken into account by means of the related sample concen- 
tration of the system c, preferable to determine experimentally and so better 
to compare with experimental results. For the sake of simplicity we assume 
the same monodisperse system as in the previous section, in order to have the 
less uncontrolled parameter-dispersity as possible. If defining c as the fraction 
of the volume occupied by all the particles (X]f^ ^0 O'^sr the total system vol- 
ume {Vt — L^), then the relationship between particle size and volume sample 
concentration is 

y ^ Vi Nv 

The orientation of the anisotropy easy axes of the particles is a parameter 
that remains fixed along the whole simulation process, and hence its value must 
be carefully selected due to its strong influence on the magnetic properties of 
the system. For studying SPM properties we have assumed a random easy-axes 
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distribution as schematized in FIGURE [71 constraint that works quite well for 
describing real systems as frozen ferrofluids or solid matrices f44l . Once the easy 
axes are placed, the simulation of the magnetic evolution of the magnetic mo- 
ments of the particles can start, whose initial orientation can also be randomly 
distributed. The magnetic properties of the particles are contained on their mag- 
netic moment pi and magnetic anisotropy under the temperature/magnetic 
field constraints applied to the system, which determine its energy. In real sys- 
tems, the evolution between an initial state with energy Ei and a final state with 
energy Ef occurs at a certain time interval, in which the event (reorientation of 
the magnetic moments in this case) can be described by a probability distribu- 
tion function. These time-dependent processes can be simulated by giving the 
system a certain amount of chances to occur, i.e. attempts to change the con- 
figuration. These attempts define the computational time and are called Monte 
Carlo steps. The MC step constitutes, therefore, the computational equivalence 
to real time units. In our simulations one MC step is defined as N trials given 
to a system of A'^-particles to change its configuration. 

The simulation of a physical process consists in varying a desired magnitude 
(temperature, magnetic field, etc) under a certain protocol and evaluating the 
energy in the new state after a certain number of MC steps, accepting or neglect- 
ing the new configuration under the chosen algorithm. For example, to simulate 
a ZFC process in a system of A^-particles the system is cooled down from a 
high temperature in regular temperature intervals every a certain fix amount of 
MC steps in zero field down to a very low temperature. In every MC step, one 
particle is selected at random and a new orientation of its magnetic moment 
is generated and accepted under the Metropolis algorithm min [l, e"^^/''^"^] . 
This is repeated TV-times in each MC step. Once the very low temperature is 
reached, a small magnetic field is applied and the process continues now while 
heating the sample [45]. The same procedure applies for simulating the mag- 
netization vs magnetic field M(H) curve, just being different the parameter to 
vary after a certain amount of MC steps. A M(H) curve is simulated by initially 
cooling the system down to the desired temperature (in zero field for our simu- 
lations), and once it is reached, a small field is applied and increased in regular 
intervals of field/MC steps up to a high field Hmax- Then, the field is decreased 
in the same manner until —Hmax is reached, and finally increased again until 
reaching once more Hmax and having completed the cycle. To illustrate the 
features of the MC steps resembling real time units, we show in Figure [TT] 
some ZFC (left panel) and M(II) (left panel) curves for fixed temperature and 
magnetic field variations, but different MC steps, corresponding to a system of 
non-interacting particles as that shown in Figure [TOl 

Left panel of Figure [TT] shows that the curves peak at lower temperatures 
for larger amount of MC steps, what reproduces the physical behavior described 
by Eq.Q: if relating the amount of MC steps with the experimental measur- 
ing time T,„, then the longer is the time interval (amount of MC steps), the 
smaller it Tg. The behavior observed in the M(H) curves of Figure [11] (right 
panel) indicates a decrease in the coercive field for larger time intervals (higher 
amount of MC steps), also as expected since the FM-like behavior represented 



17 




Figure 11: Left panel: ZFC processes simulated for a fixed temperature interval 
variation, but different MC steps (100, 200, and 500). Right panel: M(H) 
simulated processes simulated for a fixed magnetic field interval variation, but 
different MC steps (100, 150, 200, and 250). In both cases, top panel shows 
the complete processes, while bottom panel shows a magnification of the more 
remarkable aspects of the simulation (maximum of the ZFC and coercive field, 
respectively) . 

by the area in the M(H) curves is time-dependent, as discussed in Section [21 
and tends to disappear for very large times (MC steps). However, it is impor- 
tant to indicate here that, although the physical tendency coincides with the 
expected in both kind of simulations, there is not a well-established relationship 
between MC steps and real time units and some scaling must be taken in this 
aim [21] [22]. As mentioned above, time-dependent processes in systems of mag- 
netic nanoparticles are described better by means of the Landau-Lifshitz-Gilbert 
dynamic equation |19j . 

3.4 Trial computational time steps 

It turns out from the above description of the MC technique that the choosing 
of the trial new configuration with energy Ef is of primordial importance: its 
value determines the acceptation rate of the algorithm and so the velocity and 
feasibility of the simulations. In our simulations we use the so-called solid angle 
restriction scheme j27! : the new trial orientation fltriai is randomly generated 
inside a cone of angle SO around the current orientation fl. Figure [T2] illustrates 
this choice of the new trial orientation inside a cone of angle SO around the 
current orientation of the magnetic moment. [46] 
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Figure 12: Schematic drawing of the 66-cone used to generate the new trial 
orientation jltriai- 

The S9-yahie directly rules the speed of the magnetic moments' movement, 
and so its magnitude must be carefully selected in order to ensure the adequacy 
of the simulations to resemble physical processes: if 69 is too small the system 
will evolve very slowly to the quasi-equilibrium configuration, and we could be 
unable to resemble the physical process (too many MC steps would be neces- 
sary). If 59 is too large, the system can relax too fast and make the features 
we want to study to become unobservable. In Figure [TST a). we show the ZFC 
curves of the same system obtained for different values of 59. It is observed a 
tendency similar to that displayed in Figure [Til with the curves exhibiting a 
larger peak at decreasing temperature for larger (50- values. These results sup- 
port our arguments on the functioning of the trial steps as computational time 
for reproducing real time-dependent processes. 

The results plotted in Figure [13] indicate a strong dependence of the sim- 
ulated results on 59, emphasizing so that much care must be taken if trying to 
extract time- related information about the system from the simulations: those 
depend not only on the artificial MC step, but also on the trial angle chosen. 
As pointed out at the beginning of this subsection, the Monte Carlo technique 
is very successful for analyzing quasi-static processes, and so we have focused 
in this work on the more reliable aspects of the simulations, i.e., quasi-static 
analysis of complex processes with large freedom degrees. 

Figure [T3] displays also some simulations in which the d9-va,lue is con- 
sidered to be temperature-dependent, in the same way as described in Refs. 
[211 [22] . The reason for including such dependence is not related to any effort 
concerned to use realistic times from the MC steps or to reach more accurate 
time-dependence results (Figs. [TT] and [TSf a) illustrate that it constitutes a very 
complex task). Instead, our motivation was simply to provide the simulations 
with a more realistic character: it seems to us more physical to make the ther- 
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Figure 13: ZFC processes as a function of S9. (a) corresponds to different 69- 
values constant with temperature, (b) shows different temperature-dependence 
for 69, and in (c) we compare the 69=0.1=cte case, with being temperature- 
dependent as 69 = {Omty^^. In all cases, panels on the right show the temper- 
ature evolution of 69 and the 69— 0.1— cte case is shown for the comparison. 

mal fluctuations to be larger at larger temperatures. Our argument is based 
on the fact that the analogy between the time-dependence of real experiments 
and simulations is introduced by means of attempts to change the configuration, 
which are generated randomly into the cone of angle 69. Under this assumption, 
the movement of the magnetic moments is resembled by giving the system a 
certain amount of MC steps to change its configuration. 

In a real system, the magnetic moments fluctuate because of the thermal 
energy and consequently fluctuations are reduced the smaller the temperature. 
In order to reproduce this characteristic in the simulations, it appears very 
reasonable to us to include a temperature-dependence in the value of 69, since 
it is the tool used to resemble the thermal-flucuations found in real systems. 
It is worth no note here that although there is no intention of analyzing time- 
dependence in our simulations, however the temperature-dependence expression 
used has been intensively analyzed and discussed in such a context by Nowak 
et al. [21] and Cheng et al. [22]. The temperature dependence of 69 can 
be, following Refs. [211 122] and in usual computational reduced temperature 
units[3l|33l|3il3ll[35l[3i[34] t = kBT/2KV, written as 69 = (Ct)^/^, with C 
a constant value proportional to the particle inner characteristics (size, magnetic 
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moment and anisotropy) and to the dynamics of the system (gyromagnetic ratio, 
damping parameter, measuring time). The results plotted in Figure [TSTb) show 
that for values of C of the order of 0.050, not large variations are considered 
in comparison with the S0 = 0.10 shown in FIGURE [T31 We arbitrarily 

choose the 66 = 0.050t dependence value for taking into account temperature- 
dependence of 59. As a general rule, in our simulations we have considered the 
60 = 0.10 trial cone, and only very recently the temperature dependence of 69 
has been considered. For every calculation described on the next sections the 
way of choosing 69 will be properly described. 

4 Summary 

We have introduced the basic features of the so-called superparamagnetism as 
the paramagnetic-like temperature dependence that occurs in single-domain 
magnetic systems above a characteristic blocking temperature, Tb- The special 
SPM phenomena arises from the interplay between the large magnetic super- 
moment arising from the coherent rotation of the atomic magnetic moments 
in single-domain entities, and the magnetic anisotropy resultant from the crys- 
talline, shape, etc, contributions. The presence of a magnetic anisotropy defines 
a preferential orientation direction which plays a relevant energy term when the 
temperature is comparable to Tb- The direction of this anisotropy in relation 
to the applied magnetic field defines the response of the magnetic supermoment 
of the particle with temperature. In real systems, the characteristic parameters 
(mainly anisotropy and magnetic moment) vary from particle to particle due 
to the experimental difficulties to synthesize particles with perfectly controlled 
characteristics and to place them regularly (both the spatial position and the 
anisotropy orientation) . In addition, real systems are also subjected to other ef- 
fects as the interparticle magnetic dipolar interaction, adding more uncontrolled 
parameter degrees to the already per se very complex theoretical scenario. In 
this context, we introduced the Monte Carlo technique based on the Metropolis 
algorithm as a powerful tool for the study of SPM properties. 
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